This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document.
Review of the material at hand: - sumStats.df: all variables from GenTree fieldwork, averaged by population. NB plenty of missing data, together with
summary statistics derived from stairway plot analyses
See demoHistory8.Rmd for generation of values and df.
AllSp_admixtureValues.df: admixture values for K = 2 for all species.
allBio.df: bioclimatic data for all populations. See demoHistory7.Rmd.
piStats.df: all pi statistics for all populations. See demoHistory6.Rmd.
combinations of the above, more or less reduced to remove missing data.
The way forward should include at least:
filling missing data when possible
centering/standardising the variables (use scale())
merging all data sets
running full models
running models with significant variables only
recording slopes and finding a way to compare them
Start from data as recorded in ‘gtree_allspecies_env_dataset_extra_samples_20200908.csv’
bioMeanBy <- function(x,y) by(data = x, INDICES = y, FUN = mean, na.rm = T)
bioVarBy <- function(x,y) by(data = x, INDICES = y, FUN = var, na.rm = T)
library(data.table)
library(stats)
library(pspline)
library(corrr)
sourceEnvData.df <- read.csv(file = "gtree_allspecies_env_dataset_extra_samples_20200908.csv", header = T)
#these are individual-level data (mostly)
#obtaining the means:
envDataMeans.df <- data.frame(apply(X = sourceEnvData.df[c(6:9,11:99)], FUN = bioMeanBy, y = sourceEnvData.df[,5], MARGIN = 2))
#adding species names:
speciesInEnvDataMeans <- matrix(unlist(strsplit(row.names(envDataMeans.df), split = "_")),ncol = 3, byrow = T)[,2]
envDataMeans.df <-cbind(speciesInEnvDataMeans, envDataMeans.df)
names(envDataMeans.df)[1] <- "species"
#par(mfrow=c(3,3))
for (i in 2:93) hist(envDataMeans.df[,i], xlab=names(envDataMeans.df)[i], main="")
#for (i in 11:19) hist(envDataMeans.df[,i], xlab=names(envDataMeans.df)[i], main="")
##scaling the variables
#ad-hoc function to apply scale() on all data, by species
scaleBy <- function(x,y) by(data = x, INDICES = y, FUN = scale)
tmp <- apply(X = envDataMeans.df[2:94], FUN = scaleBy, y = envDataMeans.df[,1], MARGIN = 2)
## NOT working:
#envDataMeans.scaled.df <- data.frame(apply(X = envDataMeans.df[2:93], FUN = scaleBy, y = envDataMeans.df[,1], MARGIN = 2))
par(mfrow=c(2,2))
plot(x = subset(envDataMeans.df, species == "AA")$latwgs84, y = tmp$latwgs84$AA[,1])
plot(x = subset(envDataMeans.df, species == "PA")$bio11, y = tmp$bio11$PA[,1])
plot(x = subset(envDataMeans.df, species == "QP")$prec02, y = tmp$prec02$QP[,1])
plot(x = subset(envDataMeans.df, species == "PN")$tmax12, y = tmp$tmax12$PN[,1])
#everything seems to be fine
#but population names must be repositioned in front of values
#attempt one
do.call(what = c, args = tmp$latwgs84) #oh, WOW.
## AA1 AA2 AA3 AA4 AA5 AA6
## 0.660171357 0.697835260 0.588803075 0.597080470 1.556240368 1.493407739
## AA7 AA8 AA9 AA10 AA11 AA12
## -0.626350552 -0.638418636 -0.117448840 -0.866189443 -1.158301684 -0.632377955
## AA13 BP1 BP2 BP3 BP4 BP5
## -1.554451159 -0.794158017 -0.818289200 -0.087964162 0.002611020 -1.352135646
## BP6 BP7 BP8 BP9 BP10 BP11
## -1.261889468 1.773752891 1.169975198 -1.066526865 -0.984311275 -1.065563174
## BP12 BP13 BP14 BP15 BP16 BP17
## 0.218796386 0.603111484 -1.141729195 -1.328315262 0.321350737 0.268544479
## BP18 BP19 BP20 BP21 BP22 BP23
## 0.934166784 1.060550097 -0.131106538 0.178918592 0.472289480 1.567389931
## BP24 FS1 FS2 FS3 FS4 FS5
## 1.460531724 -0.169332556 -0.153105166 -0.031350933 -0.034325228 0.069042641
## FS6 FS7 FS8 FS9 FS10 FS11
## 0.067934702 0.700316055 0.694335694 -0.882002058 -0.874522106 -0.600280135
## FS12 FS13 FS14 FS15 FS16 FS17
## -0.601379044 -1.033007337 -1.032516224 0.768066118 0.785488864 -1.276609404
## FS18 FS19 FS20 FS21 FS22 FS23
## -1.278049953 -0.911476703 -0.913962592 2.197686310 2.233286637 1.592784022
## FS24 FS25 FS26 FS27 FS28 OB
## 1.677757371 -0.202872637 -0.181126692 -0.305389823 -0.305389823 NaN
## PA1 PA2 PA3 PA4 PA5 PA6
## -0.660853274 -0.660301720 -0.850016055 -0.850532441 -0.487523916 -0.487843523
## PA7 PA8 PA9 PA10 PA11 PA12
## 1.735459380 1.123146040 -1.096875801 -1.456342026 -1.462008065 -0.817863913
## PA13 PA14 PA15 PA16 PA17 PA18
## -0.821171367 0.137978376 0.317123128 0.887423939 1.097701435 1.693617100
## PA19 PA20 PA21 PA22 PA23 PA24
## -0.003792946 -0.801309371 0.213014891 0.232129054 0.052234731 1.520354632
## PA25 PC1 PC2 PC3 PC4 PC5
## 1.446251713 0.956096507 0.635036903 0.206737351 0.684857401 -0.025796084
## PC6 PC7 PC8 PC9 PC10 PH1
## 1.190225067 -1.917749758 -1.282343327 -0.630616039 0.183551977 -1.729001888
## PH2 PH3 PH4 PH5 PH6 PH7
## 0.690542909 0.738583953 -0.636895376 -0.195832516 -0.623561212 1.342337799
## PH8 PH9 PH10 PH11 PH12 PH13
## 1.229133544 -1.355144216 -0.831271434 -0.187909922 0.959418450 0.599599907
## PN1 PN2 PN3 PN4 PN5 PN6
## 1.559508463 1.560255038 -0.309592257 -0.079673711 -0.944857246 -0.948816909
## PN7 PN8 PN9 PN10 PN11 PN12
## 0.545852140 0.156108156 -1.176075978 -0.622569294 1.157427269 -0.897565672
## PO1 PO2 PO3 PO4 PO5 PO6
## -0.323857008 0.412005899 0.260625567 1.459161461 1.131590942 0.694608802
## PO7 PO8 PO9 PO10 PO11 PO12
## 0.790790880 -1.246063238 -0.895474060 0.654436277 0.497051697 -0.116804423
## PO13 PO14 PO15 PO16 PO17 PO18
## -0.130354870 -0.081381445 -0.535127969 1.003823076 1.547244964 -1.119410010
## PO19 PO20 PO21 PO22 PO23 PP1
## -0.003145183 0.002871695 -0.558607154 -0.520086326 -2.923899571 -1.185990003
## PP2 PP3 PP4 PP5 PP6 PP7
## -1.183854607 -0.212682053 -0.228728623 0.097999082 0.099465130 -0.306050336
## PP8 PP9 PP10 PP11 PP12 PP13
## -0.307437102 1.132458523 1.078708482 0.638255274 0.633746917 0.217252295
## PP14 PP15 PP16 PP17 PP18 PP19
## 0.234723108 0.595194426 0.871205102 0.860012349 0.609345214 0.611477135
## PP20 PP21 PP22 PP23 PP24 PP25
## 0.975524164 1.013528708 -2.712622232 -2.038643336 -0.344200809 -1.148686808
## PS1 PS2 PS3 PS4 PS5 PS6
## -0.707642005 -0.712297723 0.124861803 0.091191432 -1.303634631 -1.307321614
## PS7 PS8 PS9 PS10 PS11 PS12
## 1.663129573 1.104747674 -0.773561829 -0.792781313 0.634709117 0.585443967
## PS13 PS14 PS15 PS16 PS17 PS18
## -1.401168009 -1.392120280 -0.894087772 -0.893962679 0.213801175 0.386912859
## PS19 PS20 PS21 PS22 PS23 QP1
## 0.895975838 1.081569641 1.011239018 1.010935460 1.374060297 -0.270227378
## QP2 QP3 QP4 QP5 QP6 QP7
## -0.270162866 0.656147022 0.422645875 -1.135120948 -1.134960367 -0.413033126
## QP8 QP9 QP10 QP11 QP12 QP13
## -0.406099393 0.732430620 0.989767379 -0.901891398 -0.917802466 -1.683799416
## QP14 QP15 QP16 QP17 QP18 QP19
## -1.683394775 0.732428621 1.359671718 1.323005134 0.605932017 0.986667717
## QP20 TB1 TB2 TB3 TB4 TB5
## 1.007796029 0.004458195 0.085347043 -0.669433388 -0.665950293 -0.574473344
## TB6 TB7 TB8 TB9 TB10 TB11
## -0.774273977 1.012821061 -1.045168945 -0.693976538 1.831202897 1.489447290
# nesting do.call() within lapply():
tmp2 <- lapply(X = tmp, FUN = function(x) do.call(what = c, args = x))
# using do.call to fold the list into data.frame:
envDataMeans.scaled.df <- data.frame(do.call(what = cbind, args = tmp2))
#re-inserting species names in the table
speciesInEnvDataMeansScaled <- unlist(by(data = envDataMeans.df,
INDICES = envDataMeans.df$species, FUN = function(x) x$species))
envDataMeans.scaled.df <- cbind(speciesInEnvDataMeansScaled, envDataMeans.scaled.df)
names(envDataMeans.scaled.df)[1] <- "species"
#getting reordered population names from envDataMeans.df
row.names(envDataMeans.scaled.df) <- unlist(by(data = envDataMeans.df, INDICES = envDataMeans.df$species, FUN = row.names))
data.dt <- fread("all_results_stairwayplot.tsv", header = T, nrows = 1231558)
Computing smoothed curve and derivatives As a guidance, I used https://stackoverflow.com/questions/43615469/how-to-calculate-the-slope-of-a-smoothed-curve-in-r
#using birch and CH pop 5 as test data set
tmp.dt <- subset(data.dt, species == "Bpendula" & pop == "CH_BP_05")
#there are duplicated x values, keep only one of each
tmp.dt <- tmp.dt[ which(!duplicated(tmp.dt$year)),]
tmp_smoothed.curve <- smooth.Pspline(x = tmp.dt$year,
y = tmp.dt$Ne_median,
df = 20,
spar = 20,
method = 2)
tmp.f1 <- predict(tmp_smoothed.curve, x = tmp.dt$year, nderiv=1)
tmp.f2 <- predict(tmp_smoothed.curve, x = tmp.dt$year, nderiv=2)
building a generalised function
derivatives <- function(x)
{
tmp.dt <- x[which(!duplicated(x$year)),]
temp_var <- tryCatch(smooth.Pspline(x = tmp.dt$year,
y = tmp.dt$Ne_median,
df = 20,
spar = 20,
method = 2),
error=function(e) e)
ifelse(inherits(temp_var, "error"),
{
sumstats <- rep(NA, 6)
},
{
tmp_smoothed.curve <- smooth.Pspline(x = tmp.dt$year,
y = tmp.dt$theta_per_site_median,
df = 20,
spar = 20,
method = 2)
tmp.f1 <- predict(tmp_smoothed.curve, x = tmp.dt$year, nderiv=1)
tmp.f2 <- predict(tmp_smoothed.curve, x = tmp.dt$year, nderiv=2)
sumstats <- c(mean(tmp.f1), median(tmp.f1), sd(tmp.f1),mean(tmp.f2), median(tmp.f2), sd(tmp.f2))
})
names(sumstats) <- c("mean_f1", "median_f1", "sd_f1", "mean_f2", "median_f2", "sd_f2")
return(sumstats)
}
application
derivatives.sumstats <- by(data = data.dt, INDICES = data.dt$pop, FUN = derivatives, simplify = F)
derivatives.sumstats.df <- data.frame(do.call(what = rbind, args = derivatives.sumstats))
write.table(file = "derivatives_sumstats.txt", x = derivatives.sumstats.df, sep = "\t", quote = F)
adding mean, median, and sd of theta
theta_means <- by(data = data.dt, INDICES = data.dt$pop, FUN = function(x) mean(x$theta_per_site_median), simplify = F)
theta_medians <- by(data = data.dt, INDICES = data.dt$pop, FUN = function(x) median(x$theta_per_site_median), simplify = F)
theta_sd <- by(data = data.dt, INDICES = data.dt$pop, FUN = function(x) sd(x$theta_per_site_median), simplify = F)
theta_stats.df <- do.call(what = data.frame, args = list(unlist(theta_means), unlist(theta_medians), unlist(theta_sd)))
names(theta_stats.df) <- c("mean_medianTheta", "median_medianTheta", "sd_medianTheta")
all.sumstats.df <- merge(by = 0, theta_stats.df, derivatives.sumstats.df)
row.names(all.sumstats.df) <- all.sumstats.df[,1]
all.sumstats.df <- all.sumstats.df[,-1]
#adding species names:
speciesInTheta <- matrix(unlist(strsplit(row.names(all.sumstats.df), split = "_")),ncol = 3, byrow = T)[,2]
all.sumstats.df <-cbind(speciesInTheta, all.sumstats.df)
names(all.sumstats.df)[1] <- "species"
Actual scaling and centering
tmp.theta <- apply(X = all.sumstats.df[2:10], FUN = scaleBy, y = all.sumstats.df[,1], MARGIN = 2)
reassembling a table
#attempt one
# do.call(what = c, args = tmp$latwgs84) #oh, WOW.
# nesting do.call() within lapply():
tmp2.theta <- lapply(X = tmp.theta, FUN = function(x) do.call(what = c, args = x))
# using do.call to fold the list into data.frame:
all.sumstats.scaled.df <- data.frame(do.call(what = cbind, args = tmp2.theta))
#getting reordered population names from envDataMeans.df
row.names(all.sumstats.scaled.df) <- unlist(by(data = all.sumstats.df, INDICES = all.sumstats.df$species, FUN = row.names))
Importing pi stats
piStats.df <- read.table("piStats.csv", header = T)
row.names(piStats.df) <- piStats.df[,1]
Actual scaling and centering
tmp.pi <- apply(X = piStats.df[3:7], FUN = scaleBy, y = piStats.df[,2], MARGIN = 2)
reassembling a table
#attempt one
# do.call(what = c, args = tmp$latwgs84) #oh, WOW.
# nesting do.call() within lapply():
tmp2.pi <- lapply(X = tmp.pi, FUN = function(x) do.call(what = c, args = x))
# using do.call to fold the list into data.frame:
piStats.scaled.df <- data.frame(do.call(what = cbind, args = tmp2.pi))
#getting reordered population names from envDataMeans.df
row.names(piStats.scaled.df) <- unlist(by(data = piStats.df, INDICES = piStats.df$species, FUN = row.names))
NB for the time being, I focus on conStruct results with K = 2 to keep it simple NB (2) data are on External drive for the time being, so path may not be permanent
#BP
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/BP/k2/k2_conStruct.results.Robj")
BP_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/BP/k2/k2_data.block.Robj")
row.names(BP_admixtureValues.df) <- row.names(data.block$obsCov)
#FS
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/FS/k2/k2_conStruct.results.Robj")
FS_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/FS/k2/k2_data.block.Robj")
row.names(FS_admixtureValues.df) <- row.names(data.block$obsCov)
#PA
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PA/k2/k2_conStruct.results.Robj")
PA_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PA/k2/k2_data.block.Robj")
row.names(PA_admixtureValues.df) <- row.names(data.block$obsCov)
#PO
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PO/k2/k2_conStruct.results.Robj")
PO_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PO/k2/k2_data.block.Robj")
row.names(PO_admixtureValues.df) <- row.names(data.block$obsCov)
#PP
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PP/k2/k2_conStruct.results.Robj")
PP_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PP/k2/k2_data.block.Robj")
row.names(PP_admixtureValues.df) <- row.names(data.block$obsCov)
#PS
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PS/k2/k2_conStruct.results.Robj")
PS_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/PS/k2/k2_data.block.Robj")
row.names(PS_admixtureValues.df) <- row.names(data.block$obsCov)
#QP
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/QP/k2/k2_conStruct.results.Robj")
QP_admixtureValues.df <- data.frame(conStruct.results$chain_1$MAP$admix.proportions)
load("/Volumes/temporaryStorage/Documents/projects/GenTree/GenTree_T4.4/construct/QP/k2/k2_data.block.Robj")
row.names(QP_admixtureValues.df) <- row.names(data.block$obsCov)
AllSp_admixtureValues.df <- do.call(what = rbind, args = list(BP_admixtureValues.df,
FS_admixtureValues.df,
PA_admixtureValues.df,
PO_admixtureValues.df,
PP_admixtureValues.df,
PS_admixtureValues.df,
QP_admixtureValues.df))
AllSp_admixtureValues.df$PopID <- row.names(AllSp_admixtureValues.df)
#adding species names
speciesInAdmix <- matrix(unlist(strsplit(row.names(AllSp_admixtureValues.df), split = "_")),ncol = 3, byrow = T)[,2]
AllSp_admixtureValues.df <-cbind(speciesInAdmix, AllSp_admixtureValues.df)
names(AllSp_admixtureValues.df)[1] <- "species"
## NOT USED
#merging admixture info into data frames:
# View(x = merge(AllSp_admixtureValues.df, bioClimPiClean.df[,-1], all.y = T))
# bioClimPiClean.df <- merge(AllSp_admixtureValues.df, bioClimPiClean.df[,-1], all.y = T)[-c(87,88,142),]
# names(bioClimPiClean.df)[2:3] <- c("admixLayer1", "admixLayer2")
Actual scaling and centering
tmp.admix <- apply(X = AllSp_admixtureValues.df[2:3], FUN = scaleBy, y = AllSp_admixtureValues.df[,1], MARGIN = 2)
reassembling a table
#attempt one
# do.call(what = c, args = tmp$latwgs84) #oh, WOW.
# nesting do.call() within lapply():
tmp2.admix <- lapply(X = tmp.admix, FUN = function(x) do.call(what = c, args = x))
# using do.call to fold the list into data.frame:
AllSp_admixtureValues.scaled.df <- data.frame(do.call(what = cbind, args = tmp2.admix))
#getting reordered population names from AllSp_admixtureValues.df
row.names(AllSp_admixtureValues.scaled.df) <- unlist(by(data = AllSp_admixtureValues.df, INDICES = AllSp_admixtureValues.df$species, FUN = row.names))
allDivEnv.df <- merge(x = envDataMeans.scaled.df, y = AllSp_admixtureValues.scaled.df, all.y = T, by = 0)
row.names(allDivEnv.df) <- allDivEnv.df$Row.names
allDivEnv.df <- allDivEnv.df[-1]
names(allDivEnv.df)[95:96] <- c("AdmixLayer1", "AdmixLayer2")
allDivEnv.df <- merge(x = allDivEnv.df, y = piStats.scaled.df, by = 0)
row.names(allDivEnv.df) <- allDivEnv.df$Row.names
allDivEnv.df <- allDivEnv.df[-1]
allDivEnv.df <- merge(x = allDivEnv.df, y = all.sumstats.scaled.df, by = 0)
row.names(allDivEnv.df) <- allDivEnv.df$Row.names
names(allDivEnv.df)[1] <- "popId"
#performed on normalised variables
#allDivEnv.df[is.nan(allDivEnv.df)==T] <- NA
pdf(file = "normalisedDivStatsPairwisePlots.pdf", width = 16.5, height = 11.7)
plot(allDivEnv.df[98:111], pch = ".", col = "blue")
dev.off()
## quartz_off_screen
## 2
divStats.cor <- round(cor(allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98:111]),2)
corTest <- function(x,y)
{
tmp <- list()
for (i in 1:ncol(y)) tmp[[i]] <- cor.test(x,y[,i], alternative = "t", method = "p")
names(tmp) <- colnames(y)
return(tmp)
}
divStats.cortest <- apply(X = allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98:111], MARGIN = 2, FUN = corTest,
y = allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98:111])
#correlate(x = allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98:111])
# test <- corTest(x = allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98],
# y = allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98:111])
# cor.test(allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),98],
# allDivEnv.df[-c(146, which(is.na(allDivEnv.df$mean_f1))),99])
#
#
# extracting p-values
divStats.cortest_pvals <- lapply(X = divStats.cortest, FUN = function(x) lapply(X = x, function(x) x$p.value))
# identifying significant correlations with bonferroni correction
#building a matrix of significance
divStats.cortest_signif <- matrix((unlist(divStats.cortest_pvals)) < 0.05 / 14^2, ncol = 14, byrow = T)
row.names(divStats.cortest_signif) <- row.names(divStats.cor)
colnames(divStats.cortest_signif) <- colnames(divStats.cor)
#which((unlist(divStats.cortest_pvals)) < 0.05 / 14^2)
#keeping only significant values:
divStats.cortest_signif_cor <- divStats.cor * divStats.cortest_signif
divStats.cortest_signif_cor[divStats.cortest_signif_cor == 0] <- NA
write.table(x = divStats.cortest_signif_cor, file = "divStats_correlations.txt")
NB the single row for OB is removed
library(glmulti)
## Loading required package: rJava
## Loading required package: leaps
# step_piAll <- function(x)
# {
# #cat("\n********************************\n",unique(x$species),"\n")
# # init_mod <- lm(pi_all ~ ., data = x)
# # step(init_mod, scope = . ~ .^2, direction = 'forward', k = log(nrow(x)))
# glmulti(pi_all ~ latwgs84 + lonwgs84 + AdmixLayer1, data = x,
# level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
# # step(glm(formula = pi_all ~ 1, data = x),
# # scope = ~ (latwgs84 + lonwgs84 + AdmixLayer1)^2,
# # direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_piAll_T4_4", append = T)
# piAll_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_piAll)
# sink()
# #testing
# test <- glmulti(pi_all ~ latwgs84 + lonwgs84 + AdmixLayer1, data = subset(allDivEnv.df, species = "BP"),
# level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
piAll_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
piAll_by_species_T4_4.glm[[i]] <-
glmulti(pi_all ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+lonwgs84
## Crit= 64.6860050761268
## Mean crit= 70.0418366917362
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= -23.4581933274058
## Mean crit= 21.9413641543611
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+latwgs84+lonwgs84+lonwgs84:latwgs84+AdmixLayer1:lonwgs84
## Crit= 14.7683049005242
## Mean crit= 45.9578925371911
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+AdmixLayer1+lonwgs84:latwgs84+AdmixLayer1:latwgs84
## Crit= 63.0367645879122
## Mean crit= 69.0658978250577
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 41.853819703924
## Mean crit= 65.5948030215972
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+latwgs84+AdmixLayer1:latwgs84
## Crit= 59.5142191286362
## Mean crit= 67.4462304395737
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_all~1+latwgs84
## Crit= 54.413009595149
## Mean crit= 61.5923297473622
## Completed.
names(piAll_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_pi_0fold <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = pi_0fold ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_pi_0fold_T4_4", append = T)
# pi_0fold_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_pi_0fold)
# sink()
pi_0fold_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
pi_0fold_by_species_T4_4.glm[[i]] <-
glmulti(pi_0fold ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+lonwgs84
## Crit= 62.5270813109668
## Mean crit= 67.9762816233622
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= -20.0928655921797
## Mean crit= 17.5410545700456
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+latwgs84+lonwgs84+lonwgs84:latwgs84+AdmixLayer1:lonwgs84
## Crit= 18.8767969348413
## Mean crit= 49.8732904935239
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+AdmixLayer1+lonwgs84:latwgs84+AdmixLayer1:latwgs84
## Crit= 58.9901170988622
## Mean crit= 67.5515437486915
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+latwgs84+AdmixLayer1+lonwgs84:latwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 36.6167147347159
## Mean crit= 62.323815943163
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+lonwgs84+lonwgs84:latwgs84+AdmixLayer1:latwgs84
## Crit= 68.6747158993472
## Mean crit= 74.433031951216
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_0fold~1+latwgs84
## Crit= 53.9514921316646
## Mean crit= 61.2291210946449
## Completed.
names(pi_0fold_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_pi_4fold <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = pi_4fold ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_pi_4fold_T4_4", append = T)
# pi_4fold_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_pi_4fold)
# sink()
pi_4fold_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
pi_4fold_by_species_T4_4.glm[[i]] <-
glmulti(pi_4fold ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+lonwgs84
## Crit= 66.0308710003975
## Mean crit= 70.9449711458194
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= -12.8323519264042
## Mean crit= 24.592366293605
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+latwgs84+lonwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 24.0677479220453
## Mean crit= 51.9206289323742
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+AdmixLayer1+AdmixLayer1:latwgs84
## Crit= 63.1380445344821
## Mean crit= 69.3053735988875
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 48.6050061834465
## Mean crit= 67.9503592097117
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+lonwgs84+AdmixLayer1
## Crit= 68.8078518796993
## Mean crit= 74.3196311639594
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_4fold~1+latwgs84
## Crit= 54.7164677922325
## Mean crit= 61.9401379183155
## Completed.
names(pi_4fold_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_pi0_pi4 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = pi0_pi4 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_pi0_pi4_T4_4", append = T)
# pi0_pi4_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_pi0_pi4)
# sink()
pi0_pi4_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
pi0_pi4_by_species_T4_4.glm[[i]] <-
glmulti(pi0_pi4 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1+lonwgs84
## Crit= 60.5474449890675
## Mean crit= 67.946611629445
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1+lonwgs84+AdmixLayer1:latwgs84
## Crit= 41.7861376578854
## Mean crit= 49.1968348443163
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1
## Crit= 71.5081058172764
## Mean crit= 75.1341442106822
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1+AdmixLayer1+lonwgs84:latwgs84+AdmixLayer1:lonwgs84
## Crit= 26.4602338549257
## Mean crit= 50.1540257656888
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84
## Crit= 60.5627008020203
## Mean crit= 68.8911192597046
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1+latwgs84
## Crit= 68.4407798698793
## Mean crit= 75.3346170671779
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi0_pi4~1+latwgs84
## Crit= 53.3538711162489
## Mean crit= 58.6687120359456
## Completed.
names(pi0_pi4_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_pi_putNeu <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = pi_putNeu ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_pi_putNeu_T4_4", append = T)
# pi_putNeu_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_pi_putNeu)
# sink()
pi_putNeu_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
pi_putNeu_by_species_T4_4.glm[[i]] <-
glmulti(pi_putNeu ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+lonwgs84
## Crit= 64.9916673626204
## Mean crit= 70.3876099826843
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= -22.2710812480432
## Mean crit= 23.136090864477
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+latwgs84+lonwgs84+lonwgs84:latwgs84+AdmixLayer1:lonwgs84
## Crit= 14.9340119161957
## Mean crit= 44.8893078248208
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+AdmixLayer1+AdmixLayer1:latwgs84
## Crit= 63.7483895760611
## Mean crit= 69.5367257152127
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+latwgs84+AdmixLayer1+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 45.0810875045569
## Mean crit= 67.1625220816438
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+AdmixLayer1+lonwgs84:latwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 58.5212487339944
## Mean crit= 65.466846224388
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: pi_putNeu~1+latwgs84
## Crit= 54.3896209662842
## Mean crit= 61.5930032948175
## Completed.
names(pi_putNeu_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_mean_f1 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = mean_f1 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_mean_f1(theta)_T4_4", append = T)
# mean_f1_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_mean_f1)
# sink()
mean_f1_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
mean_f1_by_species_T4_4.glm[[i]] <-
glmulti(mean_f1 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1
## Crit= 61.7231399875439
## Mean crit= 68.2911056128203
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1+latwgs84
## Crit= 69.132052492742
## Mean crit= 75.2755161833566
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1
## Crit= 70.5197704201441
## Mean crit= 78.527063606074
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1+latwgs84+lonwgs84
## Crit= 48.1789687229579
## Mean crit= 55.1525997186282
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1+latwgs84+AdmixLayer1
## Crit= 51.3432509387448
## Mean crit= 55.518820089849
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1+latwgs84
## Crit= 65.5429411149059
## Mean crit= 73.0722059642032
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f1~1
## Crit= 49.918594168828
## Mean crit= 55.5847865824735
## Completed.
names(mean_f1_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_median_f1 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = median_f1 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_median_f1(theta)_T4_4", append = T)
# median_f1_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_median_f1)
# sink()
median_f1_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
median_f1_by_species_T4_4.glm[[i]] <-
glmulti(median_f1 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1
## Crit= 61.7231399875439
## Mean crit= 68.6835304760217
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1
## Crit= 73.4437265084691
## Mean crit= 79.9497816772005
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1
## Crit= 70.5197704201441
## Mean crit= 78.05420502209
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1+latwgs84+lonwgs84+AdmixLayer1:lonwgs84
## Crit= 47.5722618638525
## Mean crit= 54.2078493802
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1+latwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 50.6538477263543
## Mean crit= 55.7234998583543
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1+latwgs84
## Crit= 68.7636392421649
## Mean crit= 74.2716548010598
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f1~1+AdmixLayer1+AdmixLayer1:latwgs84
## Crit= 48.7274677120161
## Mean crit= 52.1582000861066
## Completed.
names(median_f1_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_sd_f1 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = sd_f1 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_sd_f1(theta)_T4_4", append = T)
# sd_f1_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_sd_f1)
# sink()
sd_f1_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
sd_f1_by_species_T4_4.glm[[i]] <-
glmulti(sd_f1 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1+latwgs84
## Crit= 61.6799081802282
## Mean crit= 67.5839064384881
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1
## Crit= 73.4437265084691
## Mean crit= 80.5986191095661
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1
## Crit= 70.5197704201441
## Mean crit= 75.270265993595
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1+lonwgs84+AdmixLayer1:latwgs84
## Crit= 46.0975465416833
## Mean crit= 53.8577956190487
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1+AdmixLayer1
## Crit= 49.8708841210917
## Mean crit= 55.8881003751422
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1
## Crit= 70.5197704201441
## Mean crit= 77.3042717070105
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f1~1
## Crit= 49.918594168828
## Mean crit= 55.4554257232742
## Completed.
names(sd_f1_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_mean_f2 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = mean_f2 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_mean_f2(theta)_T4_4", append = T)
# mean_f2_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_mean_f2)
# sink()
mean_f2_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
mean_f2_by_species_T4_4.glm[[i]] <-
glmulti(mean_f2 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1
## Crit= 61.7231399875439
## Mean crit= 67.9495235119085
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1+AdmixLayer1:lonwgs84
## Crit= 66.126612759353
## Mean crit= 73.122344957574
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1+AdmixLayer1:latwgs84
## Crit= 70.3165447778428
## Mean crit= 74.2267179846146
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1
## Crit= 52.8797182461919
## Mean crit= 59.9984033355274
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1+lonwgs84
## Crit= 50.3403368199975
## Mean crit= 54.8942311926022
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1+latwgs84
## Crit= 69.829666358275
## Mean crit= 75.3604761534615
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_f2~1+latwgs84+lonwgs84:latwgs84+AdmixLayer1:lonwgs84
## Crit= 46.172786800883
## Mean crit= 53.4526655301883
## Completed.
names(mean_f2_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_median_f2 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = median_f2 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_median_f2(theta)_T4_4", append = T)
# median_f2_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_median_f2)
# sink()
median_f2_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
median_f2_by_species_T4_4.glm[[i]] <-
glmulti(median_f2 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1
## Crit= 61.7231399875439
## Mean crit= 67.9773391067349
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1+latwgs84+AdmixLayer1+AdmixLayer1:lonwgs84
## Crit= 69.4541897816253
## Mean crit= 75.1662137109424
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1+lonwgs84+AdmixLayer1
## Crit= 69.9400919849175
## Mean crit= 74.8885011030612
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1+lonwgs84+AdmixLayer1:latwgs84
## Crit= 41.7262707354204
## Mean crit= 47.1530386021784
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1+latwgs84+AdmixLayer1
## Crit= 46.4784780304148
## Mean crit= 53.4392389637639
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1
## Crit= 70.5197704201441
## Mean crit= 77.7899326801508
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_f2~1+AdmixLayer1+AdmixLayer1:latwgs84
## Crit= 44.0028599589757
## Mean crit= 49.5850388071119
## Completed.
names(median_f2_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_sd_f2 <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = sd_f2 ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_sd_f2(theta)_T4_4", append = T)
# sd_f2_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_sd_f2)
# sink()
sd_f2_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
sd_f2_by_species_T4_4.glm[[i]] <-
glmulti(sd_f2 ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1
## Crit= 61.7231399875439
## Mean crit= 67.9463832692104
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1+AdmixLayer1:lonwgs84
## Crit= 66.4301940826138
## Mean crit= 73.1306722650337
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1
## Crit= 70.5197704201441
## Mean crit= 78.411378567466
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1
## Crit= 52.8797182461919
## Mean crit= 59.8668280954871
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1+latwgs84+AdmixLayer1
## Crit= 46.6686717559761
## Mean crit= 53.4719546903156
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1
## Crit= 70.5197704201441
## Mean crit= 76.3381396549914
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_f2~1+latwgs84+lonwgs84:latwgs84+AdmixLayer1:lonwgs84
## Crit= 49.3807603056406
## Mean crit= 54.1118271211957
## Completed.
names(sd_f2_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_mean_medianTheta <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = mean_medianTheta ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_mean_median(theta)_T4_4", append = T)
# mean_medianTheta_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_mean_medianTheta)
# sink()
mean_medianTheta_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
mean_medianTheta_by_species_T4_4.glm[[i]] <-
glmulti(mean_medianTheta ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1+latwgs84
## Crit= 62.8048858208063
## Mean crit= 72.3118501935249
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1+latwgs84+lonwgs84+lonwgs84:latwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 6.93857334333576
## Mean crit= 62.5926520041192
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1+lonwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 64.8882541155803
## Mean crit= 76.2364675380037
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1
## Crit= 67.5919400237546
## Mean crit= 73.1362158762125
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1+AdmixLayer1+lonwgs84:latwgs84
## Crit= 67.499988955101
## Mean crit= 76.0210168713299
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1+latwgs84+lonwgs84+AdmixLayer1:latwgs84
## Crit= 51.2770601448082
## Mean crit= 64.4133636353182
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: mean_medianTheta~1
## Crit= 58.7812650159752
## Mean crit= 65.6417923118325
## Completed.
names(mean_medianTheta_by_species_T4_4.glm) <- names(bySpecies.ls)
# step_median_medianTheta <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = median_medianTheta ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_median_median(theta)_T4_4", append = T)
# median_medianTheta_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_median_medianTheta)
# sink()
median_medianTheta_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
median_medianTheta_by_species_T4_4.glm[[i]] <-
glmulti(median_medianTheta ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1+latwgs84
## Crit= 62.1118092578636
## Mean crit= 72.0551287354967
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1+latwgs84+AdmixLayer1+AdmixLayer1:lonwgs84
## Crit= 38.8906484017752
## Mean crit= 66.0333038350318
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1+lonwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 70.7459398727705
## Mean crit= 77.8678108922215
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1
## Crit= 67.5919400237546
## Mean crit= 73.3977179380594
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1+lonwgs84
## Crit= 64.5340495752154
## Mean crit= 72.7459093604912
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1+lonwgs84
## Crit= 68.2451187691417
## Mean crit= 73.4488360662963
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: median_medianTheta~1
## Crit= 58.7812650159752
## Mean crit= 65.6585735442999
## Completed.
names(median_medianTheta_by_species_T4_4.glm) <- names(bySpecies.ls)
##NOT USED:
# step_sd_medianTheta <- function(x)
# {
# cat("\n********************************\n",unique(x$species),"\n")
# step(glm(formula = sd_medianTheta ~ 1, data = x),
# scope = ~ latwgs84 + lonwgs84 + AdmixLayer1,
# direction = "forward", k = log(nrow(x)))
# }
# sink(file = "glm_sd_median(theta)_T4_4", append = T)
# sd_medianTheta_by_species_T4_4.glm <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = step_sd_medianTheta)
# sink()
sd_medianTheta_by_species_T4_4.glm <- list()
bySpecies.ls <- by(allDivEnv.df[-146,], INDICES = factor(allDivEnv.df[-146,]$species), FUN = function(x) x)
for (i in 1:7)
{
sd_medianTheta_by_species_T4_4.glm[[i]] <-
glmulti(sd_medianTheta ~ latwgs84 + lonwgs84 + AdmixLayer1, data = bySpecies.ls[[i]],
level = 2, method = "h",fitfunction = lm, crit = 'bic', plotty = F)
}
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1
## Crit= 70.5197704201441
## Mean crit= 76.2743715643472
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1+AdmixLayer1+lonwgs84:latwgs84
## Crit= 67.3355497472131
## Mean crit= 77.045351975393
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1+latwgs84+lonwgs84+AdmixLayer1:latwgs84+AdmixLayer1:lonwgs84
## Crit= 66.1173491775956
## Mean crit= 75.9172913637316
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1
## Crit= 67.5919400237546
## Mean crit= 73.0442989157055
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1+latwgs84+lonwgs84:latwgs84
## Crit= 74.9629716485368
## Mean crit= 80.2871376954467
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1+latwgs84+lonwgs84+AdmixLayer1:latwgs84
## Crit= 54.2257643752469
## Mean crit= 64.5542011517097
## Completed.
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
##
## After 50 models:
## Best model: sd_medianTheta~1
## Crit= 58.7812650159752
## Mean crit= 65.400421015353
## Completed.
names(sd_medianTheta_by_species_T4_4.glm) <- names(bySpecies.ls)
# #where to get coefficient estimates:
# summary(lm(summary(piAll_by_species_T4_4.glm[[1]])$bestmodel,
# data=bySpecies.ls[[1]]))$coefficients[,1]
# summary(lm(summary(get(glm_T4_4[i])[[j]])$bestmodel,
# data=bySpecies.ls[[1]]))$coefficients[,1]
# collecting object names
glm_T4_4 <- ls(pattern = "_T4_4.glm")
glm_coefficients <- matrix(nrow = 14*7, ncol = 9)
k <- 0
for (i in 1:14)
{
for (j in 1:7)
{
k <- k+1
## NOT USED:
# tmp <- get(glm_T4_4[i])[[j]]$coefficients
tmp <- summary(lm(summary(get(glm_T4_4[i])[[j]])$bestmodel,
data=bySpecies.ls[[j]]))$coefficients[,1]
glm_coefficients[k,1] <- glm_T4_4[i]
glm_coefficients[k,2] <- names(get(glm_T4_4[i]))[j]
glm_coefficients[k,3] <- round(tmp[1], 3)
if(is.element("AdmixLayer1:latwgs84", names(tmp))) glm_coefficients[k,4] <- round(tmp["AdmixLayer1:latwgs84"], 3)
if(is.element("AdmixLayer1:lonwgs84", names(tmp))) glm_coefficients[k,5] <- round(tmp["AdmixLayer1:lonwgs84"], 3)
if(is.element("lonwgs84:latwgs84", names(tmp))) glm_coefficients[k,6] <- round(tmp["lonwgs84:latwgs84"], 3)
if(is.element("AdmixLayer1", names(tmp))) glm_coefficients[k,7] <- round(tmp["AdmixLayer1"], 3)
if(is.element("lonwgs84", names(tmp))) glm_coefficients[k,8] <- round(tmp["lonwgs84"], 3)
if(is.element("latwgs84", names(tmp))) glm_coefficients[k,9] <- round(tmp["latwgs84"], 3)
}
}
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
glm_coefficients.df <- data.frame(glm_coefficients)
names(glm_coefficients.df) <- c("genDivVar","species","intercept","AdmixCoef x Lat","AdmixCoef x Lon",
"Lon x Lat","AdmixCoef","Lon","Lat")
#making variable names human-friendly
glm_coefficients.df[glm_coefficients.df == "mean_f1_by_species_T4_4.glm"] <- "stairway_mean(d(theta)/dt)"
glm_coefficients.df[glm_coefficients.df == "mean_f2_by_species_T4_4.glm"] <- "stairway_mean(d2(theta)/d2t)"
glm_coefficients.df[glm_coefficients.df == "mean_medianTheta_by_species_T4_4.glm"] <- "stairway_mean(theta)"
glm_coefficients.df[glm_coefficients.df == "median_f1_by_species_T4_4.glm"] <- "stairway_median(d(theta)/dt)"
glm_coefficients.df[glm_coefficients.df == "median_f2_by_species_T4_4.glm"] <- "stairway_median(d2(theta)/d2t)"
glm_coefficients.df[glm_coefficients.df == "median_medianTheta_by_species_T4_4.glm"] <- "stairway_median(theta)"
glm_coefficients.df[glm_coefficients.df == "pi_0fold_by_species_T4_4.glm"] <- "pi_0fold"
glm_coefficients.df[glm_coefficients.df == "pi_4fold_by_species_T4_4.glm"] <- "pi_4fold"
glm_coefficients.df[glm_coefficients.df == "pi_putNeu_by_species_T4_4.glm"] <- "pi_putNeu"
glm_coefficients.df[glm_coefficients.df == "pi0_pi4_by_species_T4_4.glm"] <- "pi0/pi4"
glm_coefficients.df[glm_coefficients.df == "piAll_by_species_T4_4.glm"] <- "piAll"
glm_coefficients.df[glm_coefficients.df == "sd_f1_by_species_T4_4.glm"] <- "stairway_stdev(d(theta)/dt)"
glm_coefficients.df[glm_coefficients.df == "sd_f2_by_species_T4_4.glm"] <- "stairway_stdev(d2(theta)/d2t)"
glm_coefficients.df[glm_coefficients.df == "sd_medianTheta_by_species_T4_4.glm"] <- "stairway_stdev(theta)"
write.table(x = glm_coefficients.df, file = "glm_coefficients.txt", row.names = F, col.names = T, quote = F, sep = "\t")
# building an indicator of whether the coeeficient is significant in the full model:
# starting from glm_coefficients.df: keeping relevant lines & cols
labels.df <- glm_coefficients.df[c(1:14,36:42,64:84,92:98), c(1:2,8:9)]
# (re)creating correspondence between line names and column names in allDivEnv.df
labels.df[labels.df == "stairway_median(theta)"] <- names(allDivEnv.df)[104]
labels.df[labels.df == "stairway_mean(d(theta)/dt)"] <- names(allDivEnv.df)[106]
labels.df[labels.df == "stairway_mean(d2(theta)/d2t)"] <- names(allDivEnv.df)[109]
labels.df[labels.df == "stairway_stdev(d(theta)/dt)"] <- names(allDivEnv.df)[108]
labels.df[labels.df == "stairway_stdev(theta)"] <- names(allDivEnv.df)[105]
labels.df[labels.df == "piAll"] <- names(allDivEnv.df)[98]
labels.df[labels.df == "pi0/pi4"] <- names(allDivEnv.df)[102]
labels.df <- cbind(labels.df, glm_coefficients.df[c(1:14,36:42,64:84,92:98), 1])
names(labels.df)[5] <- "genDivVarClean"
#
pdf(file = "divStatVsGeoBySpeciesLAT.pdf", width = 20, height = 20)
par(mfrow = c(7,7))
for (i in c(98,102,104,105,106,108,109))
{
for (k in levels(as.factor(allDivEnv.df$species[-146])))
{
tmp <- subset(allDivEnv.df, species == k)[-146,]
# checking condition on significance, adding stars
tmp.labels <- subset(labels.df, species == k & genDivVar == names(tmp)[i])
if (is.na(tmp.labels[1,4]) == F) Ylab <- paste(tmp.labels[1,5], "(*)", sep = " ") else Ylab <- tmp.labels[1,5]
plot(tmp[,i] ~ tmp$latwgs84, pch = 21, bg = "grey", xaxt = "n", yaxt = "n", ylab = Ylab, xlab = "latitude", main = k,
xilm = c(min(allDivEnv.df$latwgs84), max(allDivEnv.df$latwgs84)),
ylim = c(min(allDivEnv.df[-146,i], na.rm = T), max(allDivEnv.df[-146,i], na.rm = T)))
}
}
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
dev.off()
## quartz_off_screen
## 2
#
pdf(file = "divStatVsGeoBySpeciesLON.pdf", width = 20, height = 20)
par(mfrow = c(7,7))
for (i in c(98,102,104,105,106,108,109))
{
for (k in levels(as.factor(allDivEnv.df$species[-146])))
{
tmp <- subset(allDivEnv.df, species == k)[-146,]
# checking condition on significance, adding stars
tmp.labels <- subset(labels.df, species == k & genDivVar == names(tmp)[i])
if (is.na(tmp.labels[1,3]) == F) Ylab <- paste(tmp.labels[1,5], "(*)", sep = " ") else Ylab <- tmp.labels[1,5]
plot(tmp[,i] ~ tmp$lonwgs84, pch = 21, bg = "grey", xaxt = "n", yaxt = "n", ylab = Ylab, xlab = "longitude", main = k,
xilm = c(min(allDivEnv.df$lonwgs84), max(allDivEnv.df$lonwgs84)),
ylim = c(min(allDivEnv.df[-146,i], na.rm = T), max(allDivEnv.df[-146,i], na.rm = T)))
}
}
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.window(...): "xilm" n'est pas un paramètre graphique
## Warning in plot.xy(xy, type, ...): "xilm" n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "xilm" n'est pas un
## paramètre graphique
## Warning in box(...): "xilm" n'est pas un paramètre graphique
## Warning in title(...): "xilm" n'est pas un paramètre graphique
dev.off()
## quartz_off_screen
## 2
# notice that, of course, here one cannot obtain the same slopes as in the full model used above
(using normalised values would suppress the effect)
# 1. defining functions
geoGroups <- function(x) if (is.element(x, c("PA","BP","PS"))==T) "boreal" else "temperate"
phyloGroups <- function(x) if (is.element(x, c("PA","PP","PS"))==T) "gymno" else "angio"
# reassembling relevant unscaled data
allDivGeo.raw.df <- merge(x = envDataMeans.df[1:3],piStats.df[-2], by = 0, all.y = T)
row.names(allDivGeo.raw.df)<- allDivGeo.raw.df$Row.names
allDivGeo.raw.df <- allDivGeo.raw.df[-1]
allDivGeo.raw.df <- merge(allDivGeo.raw.df, all.sumstats.df[-1], by = 0, all.y = T)
row.names(allDivGeo.raw.df)<- allDivGeo.raw.df$Row.names
allDivGeo.raw.df <- allDivGeo.raw.df[-1]
# applying functions
allDivGeo.raw.geoGroups <- sapply(X = allDivGeo.raw.df$species, FUN = geoGroups)
allDivGeo.raw.phyloGroups <- sapply(X = allDivGeo.raw.df$species, FUN = phyloGroups)
allDivGeo.raw.df <- cbind(allDivGeo.raw.geoGroups, allDivGeo.raw.phyloGroups, allDivGeo.raw.df)[-which(allDivGeo.raw.df$species == "OB"),]
names(allDivGeo.raw.df)[1:2] <- c("geoGroups","phyloGroups")
allDivGeo.raw.df$geoGroups <- factor(allDivGeo.raw.df$geoGroups)
allDivGeo.raw.df$phyloGroups <- factor(allDivGeo.raw.df$phyloGroups)
allDivGeo.raw.df$species <- factor(allDivGeo.raw.df$species)
Running models (ANOVA):
# testing
options(contrasts=c('contr.sum','contr.sum'))
contrasts(allDivGeo.raw.df$geoGroups) <- c(-0.5,0.5)
contrasts(allDivGeo.raw.df$species)
## [,1] [,2] [,3] [,4] [,5] [,6]
## BP 1 0 0 0 0 0
## FS 0 1 0 0 0 0
## PA 0 0 1 0 0 0
## PO 0 0 0 1 0 0
## PP 0 0 0 0 1 0
## PS 0 0 0 0 0 1
## QP -1 -1 -1 -1 -1 -1
summary(aov(pi_all ~ geoGroups/species, data = allDivGeo.raw.df))
## Df Sum Sq Mean Sq F value Pr(>F)
## geoGroups 1 7.000e-07 7.00e-07 18.53 2.95e-05 ***
## geoGroups:species 5 1.795e-04 3.59e-05 952.98 < 2e-16 ***
## Residuals 156 5.880e-06 4.00e-08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
effects(aov(pi_all ~ geoGroups/species, data = allDivGeo.raw.df))
## (Intercept) geoGroups1
## -3.437604e-02 -8.354193e-04
## geoGroupsboreal:species1 geoGroupstemperate:species1
## -3.199672e-03 1.134956e-02
## geoGroupstemperate:species2 geoGroupsboreal:species3
## 4.549994e-03 -4.121778e-03
## geoGroupstemperate:species4
## -1.660154e-03 -3.184523e-05
##
## 1.127230e-05 -3.845758e-04
##
## -3.815865e-04 -1.454860e-04
##
## -3.475684e-05 1.046131e-04
##
## 2.420948e-04 -1.646756e-04
##
## -1.606914e-04 -4.016890e-05
##
## -5.670911e-05 4.188323e-05
##
## 6.406339e-06 4.921999e-05
##
## 3.955862e-05 -1.122587e-04
##
## -2.562893e-04 1.500694e-05
##
## -2.338329e-05 -7.173558e-05
##
## -1.545235e-04 1.752920e-04
##
## 9.571991e-05 4.201323e-04
##
## -5.887968e-05 -4.793393e-05
##
## -2.400464e-04 -2.738791e-04
##
## 1.675410e-04 3.527013e-04
##
## 3.806537e-04 3.590371e-04
##
## 1.473222e-04 1.562395e-04
##
## 1.476028e-04 1.227030e-04
##
## 1.848759e-04 2.222841e-04
##
## -1.827190e-05 6.346906e-05
##
## 4.560770e-05 6.267649e-05
##
## -3.836787e-05 -5.606951e-07
##
## 4.310991e-04 3.624116e-04
##
## -1.133874e-04 -3.085951e-04
##
## 9.837450e-05 2.680981e-06
##
## -1.050922e-04 -5.901771e-05
##
## -5.823180e-05 5.965829e-05
##
## 7.259521e-05 -6.122726e-04
##
## 2.624607e-04 1.015907e-04
##
## 2.878589e-04 1.677725e-04
##
## 3.237854e-05 2.913715e-04
##
## 2.826790e-05 -1.320552e-05
##
## 1.539381e-05 -1.457908e-04
##
## -1.276661e-04 1.233366e-05
##
## -3.415997e-05 2.631221e-05
##
## 1.286155e-04 1.276810e-05
##
## -2.255241e-04 -3.433461e-04
##
## -3.754001e-05 -1.069760e-04
##
## -1.502020e-04 -1.308831e-04
##
## -2.739286e-04 1.688055e-04
##
## 5.125668e-05 -2.110159e-05
##
## 2.559927e-04 -1.530736e-04
##
## -1.459765e-04 1.418593e-04
##
## 1.395252e-04 -4.027643e-04
##
## -8.987361e-05 1.683631e-04
##
## 1.405098e-04 -2.720947e-04
##
## -2.337527e-04 6.986189e-05
##
## 9.121452e-05 1.964975e-04
##
## 1.421349e-04 -2.038777e-04
##
## -1.237844e-04 -1.686162e-04
##
## -1.779475e-04 -1.547538e-04
##
## -4.391958e-05 2.683465e-05
##
## -7.909757e-05 -4.911300e-04
##
## -3.167330e-04 -1.665555e-04
##
## -1.850987e-04 -6.533562e-05
##
## -1.176700e-04 1.123326e-04
##
## 1.627884e-04 -2.023571e-04
##
## -1.327043e-04 -2.945396e-04
##
## -1.781894e-04 -7.226876e-06
##
## -1.072259e-05 -9.307273e-05
##
## 4.446548e-05 -3.664236e-04
##
## -3.853351e-04 2.480223e-04
##
## 3.599447e-04 1.992351e-04
##
## -5.933843e-05 -1.419748e-04
##
## -1.060551e-05 1.171120e-04
##
## 1.336149e-04 2.218081e-04
##
## 9.906385e-05 -3.209864e-05
##
## 1.227392e-04 -1.936711e-05
##
## 2.997017e-04 2.660171e-04
##
## 4.490699e-04 3.244390e-04
##
## -3.941574e-05 -1.731951e-04
##
## 4.129331e-05 6.455946e-05
##
## -6.328115e-05 -8.475174e-05
##
## 2.939691e-04 2.901112e-04
##
## -1.100342e-04 3.717123e-04
##
## 3.967414e-07 9.266596e-06
##
## 6.018668e-05 7.527676e-05
##
## -5.216126e-05
## attr(,"assign")
## [1] 0 1 2 2 2 2 2 2 2 2 2 2 2 2
## attr(,"class")
## [1] "coef"
summary(aov(pi_all ~ phyloGroups/species, data = allDivGeo.raw.df, contrasts = "contr.sum"))
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Df Sum Sq Mean Sq F value Pr(>F)
## phyloGroups 1 1.012e-05 1.012e-05 268.6 <2e-16 ***
## phyloGroups:species 5 1.701e-04 3.402e-05 903.0 <2e-16 ***
## Residuals 156 5.880e-06 4.000e-08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(aov(pi_all ~ geoGroups/species, data = allDivGeo.raw.df))
TukeyHSD(aov(pi_all ~ geoGroups/species, data = allDivGeo.raw.df))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pi_all ~ geoGroups/species, data = allDivGeo.raw.df)
##
## $geoGroups
## diff lwr upr p adj
## temperate-boreal 0.00013197 7.140754e-05 0.0001925325 2.95e-05
##
## $`geoGroups:species`
## diff lwr upr p adj
## temperate:BP-boreal:BP NA NA NA NA
## boreal:FS-boreal:BP NA NA NA NA
## temperate:FS-boreal:BP 7.928675e-04 6.035559e-04 0.0009821790 0.0000000
## boreal:PA-boreal:BP 1.382063e-03 1.190982e-03 0.0015731430 0.0000000
## temperate:PA-boreal:BP NA NA NA NA
## boreal:PO-boreal:BP NA NA NA NA
## temperate:PO-boreal:BP -6.107229e-05 -2.582967e-04 0.0001361521 0.9986730
## boreal:PP-boreal:BP NA NA NA NA
## temperate:PP-boreal:BP -5.463787e-04 -7.374592e-04 -0.0003552982 0.0000000
## boreal:PS-boreal:BP 1.911734e-04 -3.847357e-06 0.0003861941 0.0607936
## temperate:PS-boreal:BP NA NA NA NA
## boreal:QP-boreal:BP NA NA NA NA
## temperate:QP-boreal:BP 2.999910e-03 2.794882e-03 0.0032049380 0.0000000
## boreal:FS-temperate:BP NA NA NA NA
## temperate:FS-temperate:BP NA NA NA NA
## boreal:PA-temperate:BP NA NA NA NA
## temperate:PA-temperate:BP NA NA NA NA
## boreal:PO-temperate:BP NA NA NA NA
## temperate:PO-temperate:BP NA NA NA NA
## boreal:PP-temperate:BP NA NA NA NA
## temperate:PP-temperate:BP NA NA NA NA
## boreal:PS-temperate:BP NA NA NA NA
## temperate:PS-temperate:BP NA NA NA NA
## boreal:QP-temperate:BP NA NA NA NA
## temperate:QP-temperate:BP NA NA NA NA
## temperate:FS-boreal:FS NA NA NA NA
## boreal:PA-boreal:FS NA NA NA NA
## temperate:PA-boreal:FS NA NA NA NA
## boreal:PO-boreal:FS NA NA NA NA
## temperate:PO-boreal:FS NA NA NA NA
## boreal:PP-boreal:FS NA NA NA NA
## temperate:PP-boreal:FS NA NA NA NA
## boreal:PS-boreal:FS NA NA NA NA
## temperate:PS-boreal:FS NA NA NA NA
## boreal:QP-boreal:FS NA NA NA NA
## temperate:QP-boreal:FS NA NA NA NA
## boreal:PA-temperate:FS 5.891950e-04 4.039451e-04 0.0007744450 0.0000000
## temperate:PA-temperate:FS NA NA NA NA
## boreal:PO-temperate:FS NA NA NA NA
## temperate:PO-temperate:FS -8.539398e-04 -1.045521e-03 -0.0006623589 0.0000000
## boreal:PP-temperate:FS NA NA NA NA
## temperate:PP-temperate:FS -1.339246e-03 -1.524496e-03 -0.0011539962 0.0000000
## boreal:PS-temperate:FS -6.016941e-04 -7.910057e-04 -0.0004123826 0.0000000
## temperate:PS-temperate:FS NA NA NA NA
## boreal:QP-temperate:FS NA NA NA NA
## temperate:QP-temperate:FS 2.207042e-03 2.007437e-03 0.0024066478 0.0000000
## temperate:PA-boreal:PA NA NA NA NA
## boreal:PO-boreal:PA NA NA NA NA
## temperate:PO-boreal:PA -1.443135e-03 -1.636464e-03 -0.0012498057 0.0000000
## boreal:PP-boreal:PA NA NA NA NA
## temperate:PP-boreal:PA -1.928441e-03 -2.115498e-03 -0.0017413839 0.0000000
## boreal:PS-boreal:PA -1.190889e-03 -1.381970e-03 -0.0009998087 0.0000000
## temperate:PS-boreal:PA NA NA NA NA
## boreal:QP-boreal:PA NA NA NA NA
## temperate:QP-boreal:PA 1.617847e-03 1.416563e-03 0.0018191313 0.0000000
## boreal:PO-temperate:PA NA NA NA NA
## temperate:PO-temperate:PA NA NA NA NA
## boreal:PP-temperate:PA NA NA NA NA
## temperate:PP-temperate:PA NA NA NA NA
## boreal:PS-temperate:PA NA NA NA NA
## temperate:PS-temperate:PA NA NA NA NA
## boreal:QP-temperate:PA NA NA NA NA
## temperate:QP-temperate:PA NA NA NA NA
## temperate:PO-boreal:PO NA NA NA NA
## boreal:PP-boreal:PO NA NA NA NA
## temperate:PP-boreal:PO NA NA NA NA
## boreal:PS-boreal:PO NA NA NA NA
## temperate:PS-boreal:PO NA NA NA NA
## boreal:QP-boreal:PO NA NA NA NA
## temperate:QP-boreal:PO NA NA NA NA
## boreal:PP-temperate:PO NA NA NA NA
## temperate:PP-temperate:PO -4.853064e-04 -6.786355e-04 -0.0002919773 0.0000000
## boreal:PS-temperate:PO 2.522456e-04 5.502124e-05 0.0004494701 0.0018747
## temperate:PS-temperate:PO NA NA NA NA
## boreal:QP-temperate:PO NA NA NA NA
## temperate:QP-temperate:PO 3.060982e-03 2.853857e-03 0.0032681075 0.0000000
## temperate:PP-boreal:PP NA NA NA NA
## boreal:PS-boreal:PP NA NA NA NA
## temperate:PS-boreal:PP NA NA NA NA
## boreal:QP-boreal:PP NA NA NA NA
## temperate:QP-boreal:PP NA NA NA NA
## boreal:PS-temperate:PP 7.375520e-04 5.464715e-04 0.0009286325 0.0000000
## temperate:PS-temperate:PP NA NA NA NA
## boreal:QP-temperate:PP NA NA NA NA
## temperate:QP-temperate:PP 3.546288e-03 3.345005e-03 0.0037475725 0.0000000
## temperate:PS-boreal:PS NA NA NA NA
## boreal:QP-boreal:PS NA NA NA NA
## temperate:QP-boreal:PS 2.808736e-03 2.603708e-03 0.0030137646 0.0000000
## boreal:QP-temperate:PS NA NA NA NA
## temperate:QP-temperate:PS NA NA NA NA
## temperate:QP-boreal:QP NA NA NA NA
TukeyHSD(aov(pi_all ~ phyloGroups/species, data = allDivGeo.raw.df))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pi_all ~ phyloGroups/species, data = allDivGeo.raw.df)
##
## $phyloGroups
## diff lwr upr p adj
## gymno-angio -0.0005010103 -0.0005613974 -0.0004406233 0
##
## $`phyloGroups:species`
## diff lwr upr p adj
## gymno:BP-angio:BP NA NA NA NA
## angio:FS-angio:BP 7.928675e-04 6.035559e-04 0.0009821790 0.0000000
## gymno:FS-angio:BP NA NA NA NA
## angio:PA-angio:BP NA NA NA NA
## gymno:PA-angio:BP 1.382063e-03 1.190982e-03 0.0015731430 0.0000000
## angio:PO-angio:BP -6.107229e-05 -2.582967e-04 0.0001361521 0.9986730
## gymno:PO-angio:BP NA NA NA NA
## angio:PP-angio:BP NA NA NA NA
## gymno:PP-angio:BP -5.463787e-04 -7.374592e-04 -0.0003552982 0.0000000
## angio:PS-angio:BP NA NA NA NA
## gymno:PS-angio:BP 1.911734e-04 -3.847357e-06 0.0003861941 0.0607936
## angio:QP-angio:BP 2.999910e-03 2.794882e-03 0.0032049380 0.0000000
## gymno:QP-angio:BP NA NA NA NA
## angio:FS-gymno:BP NA NA NA NA
## gymno:FS-gymno:BP NA NA NA NA
## angio:PA-gymno:BP NA NA NA NA
## gymno:PA-gymno:BP NA NA NA NA
## angio:PO-gymno:BP NA NA NA NA
## gymno:PO-gymno:BP NA NA NA NA
## angio:PP-gymno:BP NA NA NA NA
## gymno:PP-gymno:BP NA NA NA NA
## angio:PS-gymno:BP NA NA NA NA
## gymno:PS-gymno:BP NA NA NA NA
## angio:QP-gymno:BP NA NA NA NA
## gymno:QP-gymno:BP NA NA NA NA
## gymno:FS-angio:FS NA NA NA NA
## angio:PA-angio:FS NA NA NA NA
## gymno:PA-angio:FS 5.891950e-04 4.039451e-04 0.0007744450 0.0000000
## angio:PO-angio:FS -8.539398e-04 -1.045521e-03 -0.0006623589 0.0000000
## gymno:PO-angio:FS NA NA NA NA
## angio:PP-angio:FS NA NA NA NA
## gymno:PP-angio:FS -1.339246e-03 -1.524496e-03 -0.0011539962 0.0000000
## angio:PS-angio:FS NA NA NA NA
## gymno:PS-angio:FS -6.016941e-04 -7.910057e-04 -0.0004123826 0.0000000
## angio:QP-angio:FS 2.207042e-03 2.007437e-03 0.0024066478 0.0000000
## gymno:QP-angio:FS NA NA NA NA
## angio:PA-gymno:FS NA NA NA NA
## gymno:PA-gymno:FS NA NA NA NA
## angio:PO-gymno:FS NA NA NA NA
## gymno:PO-gymno:FS NA NA NA NA
## angio:PP-gymno:FS NA NA NA NA
## gymno:PP-gymno:FS NA NA NA NA
## angio:PS-gymno:FS NA NA NA NA
## gymno:PS-gymno:FS NA NA NA NA
## angio:QP-gymno:FS NA NA NA NA
## gymno:QP-gymno:FS NA NA NA NA
## gymno:PA-angio:PA NA NA NA NA
## angio:PO-angio:PA NA NA NA NA
## gymno:PO-angio:PA NA NA NA NA
## angio:PP-angio:PA NA NA NA NA
## gymno:PP-angio:PA NA NA NA NA
## angio:PS-angio:PA NA NA NA NA
## gymno:PS-angio:PA NA NA NA NA
## angio:QP-angio:PA NA NA NA NA
## gymno:QP-angio:PA NA NA NA NA
## angio:PO-gymno:PA -1.443135e-03 -1.636464e-03 -0.0012498057 0.0000000
## gymno:PO-gymno:PA NA NA NA NA
## angio:PP-gymno:PA NA NA NA NA
## gymno:PP-gymno:PA -1.928441e-03 -2.115498e-03 -0.0017413839 0.0000000
## angio:PS-gymno:PA NA NA NA NA
## gymno:PS-gymno:PA -1.190889e-03 -1.381970e-03 -0.0009998087 0.0000000
## angio:QP-gymno:PA 1.617847e-03 1.416563e-03 0.0018191313 0.0000000
## gymno:QP-gymno:PA NA NA NA NA
## gymno:PO-angio:PO NA NA NA NA
## angio:PP-angio:PO NA NA NA NA
## gymno:PP-angio:PO -4.853064e-04 -6.786355e-04 -0.0002919773 0.0000000
## angio:PS-angio:PO NA NA NA NA
## gymno:PS-angio:PO 2.522456e-04 5.502124e-05 0.0004494701 0.0018747
## angio:QP-angio:PO 3.060982e-03 2.853857e-03 0.0032681075 0.0000000
## gymno:QP-angio:PO NA NA NA NA
## angio:PP-gymno:PO NA NA NA NA
## gymno:PP-gymno:PO NA NA NA NA
## angio:PS-gymno:PO NA NA NA NA
## gymno:PS-gymno:PO NA NA NA NA
## angio:QP-gymno:PO NA NA NA NA
## gymno:QP-gymno:PO NA NA NA NA
## gymno:PP-angio:PP NA NA NA NA
## angio:PS-angio:PP NA NA NA NA
## gymno:PS-angio:PP NA NA NA NA
## angio:QP-angio:PP NA NA NA NA
## gymno:QP-angio:PP NA NA NA NA
## angio:PS-gymno:PP NA NA NA NA
## gymno:PS-gymno:PP 7.375520e-04 5.464715e-04 0.0009286325 0.0000000
## angio:QP-gymno:PP 3.546288e-03 3.345005e-03 0.0037475725 0.0000000
## gymno:QP-gymno:PP NA NA NA NA
## gymno:PS-angio:PS NA NA NA NA
## angio:QP-angio:PS NA NA NA NA
## gymno:QP-angio:PS NA NA NA NA
## angio:QP-gymno:PS 2.808736e-03 2.603708e-03 0.0030137646 0.0000000
## gymno:QP-gymno:PS NA NA NA NA
## gymno:QP-angio:QP NA NA NA NA
proportion <- table(allDivGeo.raw.df$species) / nrow(allDivGeo.raw.df)
boxplot(allDivGeo.raw.df$pi_all ~ allDivGeo.raw.df$species, col = c("blue", "green", "blue", "green", "green", "blue","green"),
width=proportion)
summary(glm(pi_all ~ latwgs84 * species, data = allDivGeo.raw.df))
##
## Call:
## glm(formula = pi_all ~ latwgs84 * species, data = allDivGeo.raw.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.705e-04 -9.182e-05 8.440e-06 9.133e-05 3.754e-04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.694e-03 1.096e-04 24.576 < 2e-16 ***
## latwgs84 1.598e-07 2.349e-06 0.068 0.94584
## species1 -9.049e-04 2.251e-04 -4.020 9.23e-05 ***
## species2 8.939e-04 2.677e-04 3.339 0.00106 **
## species3 -8.973e-04 2.143e-04 -4.186 4.83e-05 ***
## species4 -2.680e-04 3.430e-04 -0.781 0.43589
## species5 -2.815e-04 3.336e-04 -0.844 0.40016
## species6 6.313e-05 2.051e-04 0.308 0.75867
## latwgs84:species1 5.142e-06 4.362e-06 1.179 0.24034
## latwgs84:species2 -1.529e-05 5.585e-06 -2.737 0.00695 **
## latwgs84:species3 3.114e-05 4.171e-06 7.465 6.45e-12 ***
## latwgs84:species4 -9.342e-06 7.504e-06 -1.245 0.21510
## latwgs84:species5 -2.185e-05 8.011e-06 -2.727 0.00716 **
## latwgs84:species6 -9.663e-06 4.034e-06 -2.395 0.01784 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2.43043e-08)
##
## Null deviance: 1.8607e-04 on 162 degrees of freedom
## Residual deviance: 3.6213e-06 on 149 degrees of freedom
## AIC: -2379.9
##
## Number of Fisher Scoring iterations: 2
pi_all.glm <- glm(pi_all ~ latwgs84 * species, data = allDivGeo.raw.df)
# p <- ggplot(allDivGeo.raw.df, aes(x=latwgs84, y=pi_all, colour=species)) +
# geom_point(size=3) +
# #geom_line(aes(y=predict(pi_all.glm), group=Subject, size="Subjects")) +
# geom_line(data=newdat, aes(y=predict(fm2, level=0, newdata=newdat), size="Population")) +
# scale_size_manual(name="Predictions", values=c("Subjects"=0.5, "Population"=3)) +
# theme_bw(base_size=22)
#defining generic function
divGLM <- function (x, y)
{
tmp <- list()
tmp[[1]] <- list(
glm(x ~ y$latwgs84 * y$species),
lm(x[which(y$species == head(levels(y$species),1))] ~
latwgs84, data = subset(y, species == head(levels(y$species),1))),
lm(x[which(y$geoGroups == "boreal")] ~
subset(y, geoGroups == "boreal")$latwgs84),
lm(x[which(y$geoGroups == "temperate")] ~
subset(y, geoGroups == "temperate")$latwgs84),
lm(x[which(y$phyloGroups == "gymno")] ~
subset(y, phyloGroups == "gymno")$latwgs84),
lm(x[which(y$phyloGroups == "angio")] ~
subset(y, phyloGroups == "angio")$latwgs84))
names(tmp[[1]]) <- c("Sp_glm", "lastSp_lm", "geo_glm", "lastGeoGroup_glm", "phylo_glm", "lastPhyloGroup_glm")
tmp[[2]] <- list(
glm(x ~ y$lonwgs84 * y$species),
lm(x[which(y$species == head(levels(y$species),1))] ~
lonwgs84, data = subset(y, species == head(levels(y$species),1))),
lm(x[which(y$geoGroups == "boreal")] ~
subset(y, geoGroups == "boreal")$lonwgs84),
lm(x[which(y$geoGroups == "temperate")] ~
subset(y, geoGroups == "temperate")$lonwgs84),
lm(x[which(y$phyloGroups == "gymno")] ~
subset(y, phyloGroups == "gymno")$lonwgs84),
lm(x[which(y$phyloGroups == "angio")] ~
subset(y, phyloGroups == "angio")$lonwgs84))
names(tmp[[2]]) <- c("Sp_glm", "lastSp_lm", "geo_glm", "lastGeoGroup_glm", "phylo_glm", "lastPhyloGroup_glm")
return(tmp)
}
GLMs_forPlotting <- apply(X = allDivGeo.raw.df[c(7,11,13,14,15,17,18)], MARGIN = 2, FUN = divGLM, y = allDivGeo.raw.df)
Pi_all: L 965
Pi0/Pi4: L 1205
median_medianTheta: L 1450
sd_medianTheta: L 1700
mean_f1: L 1950
sd_f1: L 2200
mean_f2: L 2450
pdf(file = "Pi_all_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$pi_all,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$pi_all), max(allDivGeo.raw.df$pi_all)),
col = "blue",
xlab = "Latitude", ylab = "pi_all", main = "pi_all vs Lat")
# model for first species
curve((GLMs_forPlotting[[1]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[1]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "brown", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$pi_all,
col = "green")
curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$pi_all,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$pi_all,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$pi_all,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$pi_all,
col = "red", pch = 2)
curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$pi_all,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[1]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[1]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(57,59,65,49,33,38,35,65,60,62,35),
y = c(0.0017,0.003,0.0042,0.0015,0.002,0.0025,0.005,0.0026,0.0042,0.0033,0.003),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
# geo group models
curve((GLMs_forPlotting[[1]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[1]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[1]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[1]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_piAll_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(35,20,60,22,18,66,27,80,14,67,52),
y = c(0.002,0.003,0.0042,0.0017,0.0013,0.002,0.005,0.003,0.0027,0.0036,0.0025),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_piAll_lon.df)
pdf(file = "Pi_all_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$pi_all,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$pi_all), max(allDivGeo.raw.df$pi_all)),
col = "blue",
xlab = "Longitude", ylab = "pi_all", main = "pi_all vs Lon")
# model for first species
curve((GLMs_forPlotting[[1]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[1]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$pi_all,
col = "green")
curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$pi_all,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$pi_all,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$pi_all,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$pi_all,
col = "red", pch = 2)
curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$pi_all,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[1]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[1]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_piAll_lon.df$names,
x = labelData_piAll_lon.df$x,
y = labelData_piAll_lon.df$y,
col = labelData_piAll_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[1]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[1]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[1]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[1]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[1]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_pi0_pi4_lat.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(40,62,40,50,45,37,35,60,34,32,32),
y = c(0.2,0.24,0.35,0.275,0.3,0.45,0.23,0.31,0.26,0.31,0.24),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato")
)
View(labelData_pi0_pi4_lat.df)
pdf(file = "pi0_pi4_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$pi0_pi4,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$pi0_pi4), max(allDivGeo.raw.df$pi0_pi4)),
col = "blue",
xlab = "Latitude", ylab = "pi0_pi4", main = "pi0_pi4 vs Lat")
# model for first species
curve((GLMs_forPlotting[[2]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[2]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$pi0_pi4,
col = "green")
curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$pi0_pi4,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$pi0_pi4,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$pi0_pi4,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$pi0_pi4,
col = "red", pch = 2)
curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$pi0_pi4,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[2]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[2]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_pi0_pi4_lat.df$names,
x = labelData_pi0_pi4_lat.df$x,
y = labelData_pi0_pi4_lat.df$y,
col = labelData_pi0_pi4_lat.df$col)
# geo group models
curve((GLMs_forPlotting[[2]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[2]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[2]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[2]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_pi0_pi4_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(35,-5,62,-8,18,50,27,80,24,50,40),
y = c(0.19,0.225,0.35,0.26,0.29,0.44,0.24,0.32,0.23,0.38,0.22),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_pi0_pi4_lon.df)
pdf(file = "pi0_pi4_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$pi0_pi4,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$pi0_pi4), max(allDivGeo.raw.df$pi0_pi4)),
col = "blue",
xlab = "Longitude", ylab = "pi0_pi4", main = "pi0_pi4 vs Lon")
# model for first species
curve((GLMs_forPlotting[[2]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[2]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$pi0_pi4,
col = "green")
curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$pi0_pi4,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$pi0_pi4,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$pi0_pi4,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$pi0_pi4,
col = "red", pch = 2)
curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$pi0_pi4,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[2]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[2]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_pi0_pi4_lon.df$names,
x = labelData_pi0_pi4_lon.df$x,
y = labelData_pi0_pi4_lon.df$y,
col = labelData_pi0_pi4_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[2]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[2]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[2]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[2]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[2]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_median_medianTheta_lat.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(67,50,62,50,40,66,36,38,33,32,32),
y = c(0.007,0.0035,0.0085,0.001,0.0012,0.0028,0.009,0.005,0.001,0.0025,0.0035),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato")
)
View(labelData_median_medianTheta_lat.df)
pdf(file = "median_medianTheta_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$median_medianTheta,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$median_medianTheta), max(allDivGeo.raw.df$median_medianTheta)),
col = "blue",
xlab = "Latitude", ylab = "median_medianTheta", main = "median_medianTheta vs Lat")
# model for first species
curve((GLMs_forPlotting[[3]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[3]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$median_medianTheta,
col = "green")
curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$median_medianTheta,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$median_medianTheta,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$median_medianTheta,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$median_medianTheta,
col = "red", pch = 2)
curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$median_medianTheta,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[3]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[3]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_median_medianTheta_lat.df$names,
x = labelData_median_medianTheta_lat.df$x,
y = labelData_median_medianTheta_lat.df$y,
col = labelData_median_medianTheta_lat.df$col)
# geo group models
curve((GLMs_forPlotting[[3]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[3]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[3]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[3]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_median_medianTheta_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(75,19,40,15,5,50,0,80,23,70,-8),
y = c(0.0045,0.003,0.0085,0.002,0.0012,0.0015,0.0095,0.0055,0.0042,0.0075,0.0047),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_median_medianTheta_lon.df)
pdf(file = "median_medianTheta_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$median_medianTheta,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$median_medianTheta), max(allDivGeo.raw.df$median_medianTheta)),
col = "blue",
xlab = "Longitude", ylab = "median_medianTheta", main = "median_medianTheta vs Lon")
# model for first species
curve((GLMs_forPlotting[[3]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[3]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$median_medianTheta,
col = "green")
curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$median_medianTheta,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$median_medianTheta,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$median_medianTheta,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$median_medianTheta,
col = "red", pch = 2)
curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$median_medianTheta,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[3]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[3]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_median_medianTheta_lon.df$names,
x = labelData_median_medianTheta_lon.df$x,
y = labelData_median_medianTheta_lon.df$y,
col = labelData_median_medianTheta_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[3]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[3]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[3]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[3]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[3]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_sd_medianTheta_lat.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(0,0,0,0,0,38,0,38,0,38,0),
y = c(0,0,0,0,0,0.0175,0,0.0075,0,0.006,0),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato")
)
View(labelData_sd_medianTheta_lat.df)
pdf(file = "sd_medianTheta_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$sd_medianTheta,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$sd_medianTheta), max(allDivGeo.raw.df$sd_medianTheta)),
col = "blue",
xlab = "Latitude", ylab = "sd_medianTheta", main = "sd_medianTheta vs Lat")
# model for first species
curve((GLMs_forPlotting[[4]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[4]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$sd_medianTheta,
col = "green")
curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$sd_medianTheta,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$sd_medianTheta,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$sd_medianTheta,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$sd_medianTheta,
col = "red", pch = 2)
curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$sd_medianTheta,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[4]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[4]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_sd_medianTheta_lat.df$names,
x = labelData_sd_medianTheta_lat.df$x,
y = labelData_sd_medianTheta_lat.df$y,
col = labelData_sd_medianTheta_lat.df$col)
# geo group models
curve((GLMs_forPlotting[[4]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[4]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[4]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[4]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_sd_medianTheta_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(100,100,100,100,100,38,100,0,100,0,100),
y = c(0,0,0,0,0,0.0125,0,0.009,0,0.005,0),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_sd_medianTheta_lon.df)
pdf(file = "sd_medianTheta_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$sd_medianTheta,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$sd_medianTheta), max(allDivGeo.raw.df$sd_medianTheta)),
col = "blue",
xlab = "Longitude", ylab = "sd_medianTheta", main = "sd_medianTheta vs Lon")
# model for first species
curve((GLMs_forPlotting[[4]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[4]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$sd_medianTheta,
col = "green")
curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$sd_medianTheta,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$sd_medianTheta,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$sd_medianTheta,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$sd_medianTheta,
col = "red", pch = 2)
curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$sd_medianTheta,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[4]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[4]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_sd_medianTheta_lon.df$names,
x = labelData_sd_medianTheta_lon.df$x,
y = labelData_sd_medianTheta_lon.df$y,
col = labelData_sd_medianTheta_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[4]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[4]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[4]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[4]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[4]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_mean_f1_lat.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(65,40,40,32,32,65,36,50,33,50,33),
y = c(1e-8,3e-8,-0.5e-8,0,-2.5e-8,-6e-8,-0.25e-8,-1.75e-8,-1.25e-8,-3.75e-8,2e-8),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato")
)
View(labelData_mean_f1_lat.df)
pdf(file = "mean_f1_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$mean_f1,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$mean_f1, na.rm = T), max(allDivGeo.raw.df$mean_f1, na.rm = T)),
col = "blue",
xlab = "Latitude", ylab = "mean_d(theta)/dt", main = "mean_d(theta)/dt vs Lat")
# model for first species
curve((GLMs_forPlotting[[5]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[5]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$mean_f1,
col = "green")
curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$mean_f1,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$mean_f1,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$mean_f1,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$mean_f1,
col = "red", pch = 2)
curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$mean_f1,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[5]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[5]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_mean_f1_lat.df$names,
x = labelData_mean_f1_lat.df$x,
y = labelData_mean_f1_lat.df$y,
col = labelData_mean_f1_lat.df$col)
# geo group models
curve((GLMs_forPlotting[[5]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[5]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[5]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[5]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_mean_f1_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(80,20,50,-5,5,65,-7,80,33,40,40),
y = c(1e-8,3e-8,-0.75e-8,1e-8,-4e-8,-6e-8,-0.25e-8,-1.5e-8,1.25e-8,-2.75e-8,0.75e-8),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_mean_f1_lon.df)
pdf(file = "mean_f1_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$mean_f1,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$mean_f1, na.rm = T), max(allDivGeo.raw.df$mean_f1, na.rm = T)),
col = "blue",
xlab = "Longitude", ylab = "mean_d(theta)/dt", main = "mean_d(theta)/dt vs Lon")
# model for first species
curve((GLMs_forPlotting[[5]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[5]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$mean_f1,
col = "green")
curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$mean_f1,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$mean_f1,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$mean_f1,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$mean_f1,
col = "red", pch = 2)
curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$mean_f1,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[5]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[5]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_mean_f1_lon.df$names,
x = labelData_mean_f1_lon.df$x,
y = labelData_mean_f1_lon.df$y,
col = labelData_mean_f1_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[5]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[5]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[5]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[5]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[5]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_sd_f1_lat.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(0,60,0,32,37,50,36,62,35,50,33),
y = c(1e-8,5e-8,-0.5e-8,0,6.5e-8,8e-8,-0.25e-8,2.75e-8,3.5e-8,5e-8,2e-8),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato")
)
View(labelData_sd_f1_lat.df)
pdf(file = "sd_f1_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$sd_f1,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$sd_f1, na.rm = T), max(allDivGeo.raw.df$sd_f1, na.rm = T)),
col = "blue",
xlab = "Latitude", ylab = "sd_d(theta)/dt", main = "sd_d(theta)/dt vs Lat")
# model for first species
curve((GLMs_forPlotting[[6]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[6]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$sd_f1,
col = "green")
curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$sd_f1,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$sd_f1,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$sd_f1,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$sd_f1,
col = "red", pch = 2)
curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$sd_f1,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[6]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[6]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_sd_f1_lat.df$names,
x = labelData_sd_f1_lat.df$x,
y = labelData_sd_f1_lat.df$y,
col = labelData_sd_f1_lat.df$col)
# geo group models
curve((GLMs_forPlotting[[6]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[6]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[6]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[6]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_sd_f1_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(50,5,30,-10,-10,50,15,80,30,50,70),
y = c(0,4.5e-8,1e-8,1e-8,4e-8,8e-8,-0.25e-8,2.75e-8,2e-8,4e-8,1.5e-8),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_sd_f1_lon.df)
pdf(file = "sd_f1_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$sd_f1,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$sd_f1, na.rm = T), max(allDivGeo.raw.df$sd_f1, na.rm = T)),
col = "blue",
xlab = "Longitude", ylab = "sd_d(theta)/dt", main = "sd_d(theta)/dt vs Lon")
# model for first species
curve((GLMs_forPlotting[[6]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[6]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$sd_f1,
col = "green")
curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$sd_f1,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$sd_f1,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$sd_f1,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$sd_f1,
col = "red", pch = 2)
curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$sd_f1,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[6]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[6]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_sd_f1_lon.df$names,
x = labelData_sd_f1_lon.df$x,
y = labelData_sd_f1_lon.df$y,
col = labelData_sd_f1_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[6]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[6]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[6]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[6]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[6]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_mean_f2_lat.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(0,40,0,32,37,50,36,62,35,50,33),
y = c(1e-8,-0.5e-10,-0.5e-8,0,6.5e-8,8e-8,-0.25e-8,2.75e-8,3.5e-8,5e-8,2e-8),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato")
)
View(labelData_mean_f2_lat.df)
pdf(file = "mean_f2_latBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$mean_f2,
xlim = c(min(allDivGeo.raw.df$latwgs84), max(allDivGeo.raw.df$latwgs84)),
ylim = c(min(allDivGeo.raw.df$mean_f2, na.rm = T), max(allDivGeo.raw.df$mean_f2, na.rm = T)),
col = "blue",
xlab = "Latitude", ylab = "mean_d2(theta)/d2t", main = "mean_d2(theta)/d2t vs Lat")
# model for first species
curve((GLMs_forPlotting[[7]][[1]][[2]]$coefficients["latwgs84"])*x+
GLMs_forPlotting[[7]][[1]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$latwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$mean_f2,
col = "green")
curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$latwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$mean_f2,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$latwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$mean_f2,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$latwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$mean_f2,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$latwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$latwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$mean_f2,
col = "red", pch = 2)
curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$latwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$latwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$mean_f2,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$latwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[7]][[1]][[1]]$coefficients["y$latwgs84"])*x+
# GLMs_forPlotting[[7]][[1]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$latwgs84),
# to = max(allDivGeo.raw.df$latwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_mean_f2_lat.df$names,
x = labelData_mean_f2_lat.df$x,
y = labelData_mean_f2_lat.df$y,
col = labelData_mean_f2_lat.df$col)
# geo group models
curve((GLMs_forPlotting[[7]][[1]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[1]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$latwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[7]][[1]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[1]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$latwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[7]][[1]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[1]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$latwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[7]][[1]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[1]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$latwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
labelData_mean_f2_lon.df <- data.frame(
names = c(levels(allDivGeo.raw.df$species),levels(allDivGeo.raw.df$geoGroups),levels(allDivGeo.raw.df$phyloGroups)[2:1]),
x = c(100,0,30,-10,-10,50,15,80,30,50,70),
y = c(0,-1e-10,1e-8,1e-8,4e-8,8e-8,-0.25e-8,2.75e-8,2e-8,4e-8,1.5e-8),
col = c("blue", "green", "darkgreen", "orange", "purple", "red", "brown","navyblue","navyblue","tomato","tomato"))
View(labelData_mean_f2_lon.df)
pdf(file = "mean_f2_lonBYspeciesAndGeoGroup.pdf")
plot(x = subset(allDivGeo.raw.df, species == "BP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "BP")$mean_f2,
xlim = c(min(allDivGeo.raw.df$lonwgs84), max(allDivGeo.raw.df$lonwgs84)),
ylim = c(min(allDivGeo.raw.df$mean_f2, na.rm = T), max(allDivGeo.raw.df$mean_f2, na.rm = T)),
col = "blue",
xlab = "Longitude", ylab = "mean_d2(theta)/d2t", main = "mean_d2(theta)/d2t vs Lon")
# model for first species
curve((GLMs_forPlotting[[7]][[2]][[2]]$coefficients["lonwgs84"])*x+
GLMs_forPlotting[[7]][[2]][[2]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "BP")$lonwgs84),
lty = "dashed", col = "blue", add = T)
# dots and modesl for all other species
points(x = subset(allDivGeo.raw.df, species == "FS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "FS")$mean_f2,
col = "green")
curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesFS"])*x+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$speciesFS"],
from = min(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "FS")$lonwgs84),
lty = "dashed", col = "green", add = T)
points(x = subset(allDivGeo.raw.df, species == "PA")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PA")$mean_f2,
col = "darkgreen", pch = 0)
curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPA"])*x+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$speciesPA"],
from = min(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PA")$lonwgs84),
lty = "dashed", col = "darkgreen", add = T)
points(x = subset(allDivGeo.raw.df, species == "PO")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PO")$mean_f2,
col = "orange", pch = 0)
curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPO"])*x+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$speciesPO"],
from = min(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PO")$lonwgs84),
lty = "dashed", col = "orange", add = T)
points(x = subset(allDivGeo.raw.df, species == "PP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PP")$mean_f2,
col = "purple", pch = 2)
curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPP"])*x+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$speciesPP"],
from = min(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PP")$lonwgs84),
lty = "dashed", col = "purple", add = T)
points(x = subset(allDivGeo.raw.df, species == "PS")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "PS")$mean_f2,
col = "red", pch = 2)
curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesPS"])*x+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$speciesPS"],
from = min(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "PS")$lonwgs84),
lty = "dashed", col = "red", add = T)
points(x = subset(allDivGeo.raw.df, species == "QP")$lonwgs84,
y = subset(allDivGeo.raw.df, species == "QP")$mean_f2,
col = "brown", pch = 0)
curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84:y$speciesQP"])*x+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"]+
GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$speciesQP"],
from = min(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
to = max(subset(allDivGeo.raw.df, species == "QP")$lonwgs84),
lty = "dashed", col = "brown", add = T)
## NOT USED
# overall mean
# curve((GLMs_forPlotting[[7]][[2]][[1]]$coefficients["y$lonwgs84"])*x+
# GLMs_forPlotting[[7]][[2]][[1]]$coefficients["(Intercept)"],
# from = min(allDivGeo.raw.df$lonwgs84),
# to = max(allDivGeo.raw.df$lonwgs84),
# lwd = 2, col = "lightgrey", add = T)
text(labels = labelData_mean_f2_lon.df$names,
x = labelData_mean_f2_lon.df$x,
y = labelData_mean_f2_lon.df$y,
col = labelData_mean_f2_lon.df$col)
# geo group models
curve((GLMs_forPlotting[[7]][[2]][[3]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[2]][[3]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "boreal")$lonwgs84),
col = "navyblue", add = T)
curve((GLMs_forPlotting[[7]][[2]][[4]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[2]][[4]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
to = max(subset(allDivGeo.raw.df, geoGroups == "temperate")$lonwgs84),
col = "navyblue", lty = "dotted", add = T)
# phylogroup models
curve((GLMs_forPlotting[[7]][[2]][[5]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[2]][[5]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "gymno")$lonwgs84),
col = "tomato", add = T)
curve((GLMs_forPlotting[[7]][[2]][[6]]$coefficients[2])*x+
GLMs_forPlotting[[7]][[2]][[6]]$coefficients["(Intercept)"],
from = min(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
to = max(subset(allDivGeo.raw.df, phyloGroups == "angio")$lonwgs84),
col = "tomato", lty = "dotted", add = T)
dev.off()
## quartz_off_screen
## 2
Notice that this also needs to be carried out on raw estimates, so as to be able to detect differences in mean between groups, in addition to differences in slope
# defining functions
geoBYgeoGroup_GLM <- function(x)
{
tmp1 <- glm(x ~ allDivGeo.raw.df$latwgs84 * allDivGeo.raw.df$geoGroups)
tmp2 <- glm(x ~ allDivGeo.raw.df$lonwgs84 * allDivGeo.raw.df$geoGroups)
tmp3 <- TukeyHSD(x = aov(glm(x ~ allDivGeo.raw.df$latwgs84 * allDivGeo.raw.df$geoGroups)), which = "allDivGeo.raw.df$geoGroups", suppressWarnings = T)
tmp4 <- TukeyHSD(x = aov(glm(x ~ allDivGeo.raw.df$lonwgs84 * allDivGeo.raw.df$geoGroups)), which = "allDivGeo.raw.df$geoGroups")
return(list(tmp1,tmp2,tmp3,tmp4))
}
geoBYphyloGroup_GLM <- function(x)
{
tmp1 <- glm(x ~ allDivGeo.raw.df$latwgs84 * allDivGeo.raw.df$phyloGroups)
tmp2 <- glm(x ~ allDivGeo.raw.df$lonwgs84 * allDivGeo.raw.df$phyloGroups)
tmp3 <- TukeyHSD(x = aov(glm(x ~ allDivGeo.raw.df$latwgs84 * allDivGeo.raw.df$phyloGroups)), which = "allDivGeo.raw.df$phyloGroups")
tmp4 <- TukeyHSD(x = aov(glm(x ~ allDivGeo.raw.df$lonwgs84 * allDivGeo.raw.df$phyloGroups)), which = "allDivGeo.raw.df$phyloGroups")
return(list(tmp1,tmp2,tmp3,tmp4))
}
#running functions
geoBYgeoGroup.glm <- apply(X = allDivGeo.raw.df[c(7,11,13,14,15,17,18)], MARGIN = 2, FUN = geoBYgeoGroup_GLM)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$geoGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$geoGroups
geoBYphyloGroup.glm <- apply(X = allDivGeo.raw.df[c(7,11,13,14,15,17,18)], MARGIN = 2, FUN = geoBYphyloGroup_GLM)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$latwgs84, allDivGeo.raw.df$phyloGroups
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## allDivGeo.raw.df$lonwgs84, allDivGeo.raw.df$phyloGroups
# threshold for alpha, Bonferroni-corrected
thr <- 0.05/28
# setting up matrix
geoBySpGroup.mt <- matrix(nrow = 7, ncol = 9)
colnames(geoBySpGroup.mt) <- c("divStat",
"LATgeoGroupDiff","LATgeoGroupInteraction",
"LONgeoGroupDiff","LONgeoGroupInteraction",
"LATphyloGroupDiff","LATphyloGroupInteraction",
"LONphyloGroupDiff","LONphyloGroupInteraction")
# collecting data, conditional on significance
for (i in 1:7)
{
geoBySpGroup.mt[i,1] <- names(geoBYgeoGroup.glm)[i]
if (geoBYgeoGroup.glm[[i]][[3]][[1]][4] < thr) geoBySpGroup.mt[i,2] <- -geoBYgeoGroup.glm[[i]][[3]][[1]][1]
if (summary(geoBYgeoGroup.glm[[i]][[1]])$coefficients[4,4] < thr) geoBySpGroup.mt[i,3] <- -geoBYgeoGroup.glm[[i]][[1]]$effects[4]
if (geoBYgeoGroup.glm[[i]][[4]][[1]][4] < thr) geoBySpGroup.mt[i,4] <- -geoBYgeoGroup.glm[[i]][[4]][[1]][1]
if (summary(geoBYgeoGroup.glm[[i]][[2]])$coefficients[4,4] < thr) geoBySpGroup.mt[i,5] <- -geoBYgeoGroup.glm[[i]][[2]]$effects[4]
if (geoBYphyloGroup.glm[[i]][[3]][[1]][4] < thr) geoBySpGroup.mt[i,6] <- geoBYphyloGroup.glm[[i]][[3]][[1]][1]
if (summary(geoBYphyloGroup.glm[[i]][[1]])$coefficients[4,4] < thr) geoBySpGroup.mt[i,7] <- geoBYphyloGroup.glm[[i]][[1]]$effects[4]
if (geoBYphyloGroup.glm[[i]][[4]][[1]][4] < thr) geoBySpGroup.mt[i,8] <- geoBYphyloGroup.glm[[i]][[4]][[1]][1]
if (summary(geoBYphyloGroup.glm[[i]][[2]])$coefficients[4,4] < thr) geoBySpGroup.mt[i,9] <- geoBYphyloGroup.glm[[i]][[2]]$effects[4]
}
colnames(geoBySpGroup.mt)[c(2,4,6,8)] <- c("mean(boreal)-mean(temp)_LAT","mean(boreal)-mean(temp)_LON","mean(gymno)-mean(angio)_LAT","mean(gymno)-mean(angio)_LON")
colnames(geoBySpGroup.mt)[c(3,5,7,9)] <- c("slope(boreal)-slope(temp)_LAT","slope(boreal)-slope(temp)_LON","slope(gymno)-slope(angio)_LAT","slope(gymno)-slope(angio)_LON")
geoBySpGroup.mt[,2:9] <- signif(as.numeric(geoBySpGroup.mt[,2:9]),3)
write.table(x = geoBySpGroup.mt, file = "geoBySpGroup.txt", quote = F)